R包vegan的群落对应分析（CA）

对应分析（Correspondence Analysis，CA）在群落分析中常用于物种组成数据，排序结果展示的是样方间χ²距离。χ²距离不受双零问题的影响，适用于原始物种多度数据（和PCA相比，物种多度数据无需转化）。

CA分析对象必须是频度或类频度、同量纲的非负数据，如物种个体计数或有-无数据。

对于CA的基本概念描述及其在群落分析中的应用，可参考前文。

本篇以某16S扩增子测序所得细菌群落数据，简介在R中执行CA的过程。

示例文件、R脚本等，已上传至百度盘：

https://pan.baidu.com/s/1YEsOIxgk8VWEjJ5KDRdLBA

文件“phylum_table.txt”为通过16S测序所得的细菌类群丰度表格，统计至细菌“门水平”（即界门纲目科属种的“门”这一水平）。我未使用OTU水平的，仅方便用于CA的展示；并且，作为示例，暂且忽略本数据更适用线性模型还是单峰模型

CA分析

通过CA，查看数据中各样方的物种组成差异，或者物种分布状态。

CA执行

vegan包中，可通过cca()执行CA，该函数的输入格式大致如下。

library(vegan)

#详情 ?cca
#带协变量的 CCA
cca(Y, X, W)
#不带协变量的 CCA，即常规 CCA
cca(Y, X)
#CA
cca(Y)

其中，Y为响应变量矩阵，X为解释变量矩阵，W为偏CCA分析需要输入的协变量矩阵。X或W是在CCA分析中才需要提供的，对于CA而言，只需提供响应变量矩阵Y即可，即由对象（行）和变量（列）组成的矩阵。

这里Y就是物种组成矩阵了，方法如下（和PCA的函数用法很相似）。

#读取物种数据
phylum <- read.delim('phylum_table.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)

##CA 排序，详情 ?cca
ca_phylum <- cca(phylum)

#I 型标尺
summary(ca_phylum, scaling = 1)
#II 型标尺
#summary(ca_phylum, scaling = 2)

CA结果解读

类似于PCA，CA也有两种标尺选择，通过summary()展示结果时指定参数“scaling=1或2”。I型标尺图中样方的点是物种多度的形心，更适合解释样方之间的关系和样方的梯度排列。II型标尺图中物种的点是样方的形心，更适合解释物种之间的关系和梯度分布。

关于对CA轴的特征，以及I型和II型标尺的基本描述，可参考前文。

这里展示了按I型标尺缩放后排序结果（参数scaling = 1），部分内容如下所示。

主要内容提取

可以将重要的信息提取出来，例如：

#各 CA 轴的特征根，特征根是固定的，和标尺选择无关
ca_eig <- ca_phylum$CA$eig
#除以特征根总和，即可得各 CA 轴的解释量
ca_exp <- ca_phylum$CA$eig / sum(ca_phylum$CA$eig)

#提取样方得分，以 I 型标尺为例，以前两轴为例
#scores() 提取样方坐标
site.scaling1 <- scores(ca_phylum, choices = 1:2, scaling = 1, display = 'site')
#或者在 summary 中提取
site.scaling1 <- summary(ca_phylum, scaling = 1)$sites[ ,1:2]
site.scaling1
#write.table(site.scaling1, 'site.scaling1.txt', col.names = NA, sep = '\t', quote = FALSE)

#提取物种得分，以 II 型标尺为例，以前两轴为例
#scores() 提取物种坐标
phylum.scaling2 <- scores(ca_phylum, choices = 1:2, scaling = 2, display = 'sp')
#或者 summary 中提取
phylum.scaling2 <- summary(ca_phylum, scaling = 2)$species[ ,1:2]
phylum.scaling2
#write.table(phylum.scaling2, 'phylum.scaling2.txt', col.names = NA, sep = '\t', quote = FALSE)

CA排序图

最后根据CA所得的样方或物种的排序坐标，绘图就可以了。CA排序图的解读方式，可参考前文。

通常，样方和物种直接在对应坐标处绘制为点。（与PCA的展示略有不同，PCA一般将变量展示为箭头向量形式，CA一般直接将变量以点表示）

#ordiplot() 简洁版作图观测，详情 ?ordiplot
#I 型标尺，样方点是物种点的形心
#II 型标尺，物种点是样方点的形心
ca1 <- paste('CA1:', round(ca_exp[1]*100, 2), '%')
ca2 <- paste('CA2:',round(ca_exp[2]*100, 2), '%')

par(mfrow = c(2, 2))
ordiplot(ca_phylum, display = 'site', type = 'text', choices = c(1, 2), scaling = 1, main = 'I 型标尺，仅样方', xlab = ca1, ylab = ca2)
ordiplot(ca_phylum, display = 'sp', type = 'text', choices = c(1, 2), scaling = 2, main = 'II 型标尺，仅物种', xlab = ca1, ylab = ca2)
ordiplot(ca_phylum, type = 'point', choices = c(1, 2), scaling = 1, main = 'I 型标尺，双序图', xlab = ca1, ylab = ca2)
ordiplot(ca_phylum, type = 'point', choices = c(1, 2), scaling = 2, main = 'II 型标尺，双序图', xlab = ca1, ylab = ca2)
#ordiplot() 作图，可视化调整后，这里以 I 型标尺为例，详情 ?ordiplot
#仅样方
png('ca_phylum.png', width = 600, height = 600, res = 100, units = 'px')
ordiplot(ca_phylum, type = 'n', display = 'site', choices = c(1, 2), scaling = 1, main = 'I 型标尺，仅样方', xlab = ca1, ylab = ca2)
points(ca_phylum, display = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 4), rep('orange', 4), rep('green3', 4)), col = NA, cex = 1.2)
dev.off()

#带物种
#ordiplot(ca_phylum, type = 'n', choices = c(1, 2), scaling = 1, main = 'I 型标尺，双序图', xlab = ca1, ylab = ca2)
#points(ca_phylum, display = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 4), rep('orange', 4), rep('green3', 4)), col = NA, cex = 1.2)
#text(ca_phylum, display = 'sp', choices = c(1, 2), scaling = 1, col = 'blue', cex = 0.8)
#ggplot2 作图，I 型标尺为例，展示所有样方以及部分物种（top10 丰度门）
#提取样方和物种排序坐标，前两轴
site.scaling1 <- data.frame(summary(ca_phylum, scaling = 1)$sites[ ,1:2])
phylum.scaling1 <- data.frame(summary(ca_phylum, scaling = 1)$species[ ,1:2])

#添加分组信息
site.scaling1$sample <- rownames(site.scaling1)
site.scaling1$group <- c(rep('A', 4), rep('B', 4), rep('C', 4))
phylum.scaling1$group <- rownames(phylum.scaling1)

#只保留 top10 丰度的细菌门
top10_phylum <- names(colSums(phylum)[order(colSums(phylum), decreasing = TRUE)][1:10])
top10_phylum.scaling1 <- phylum.scaling1[top10_phylum, ]

#ggplot2 作图
library(ggplot2)
library(ggrepel) #用于 geom_text_repel() 添加标签

p <- ggplot(site.scaling1, aes(CA1, CA2)) +
geom_point(aes(color = group)) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent'), plot.title = element_text(hjust = 0.5)) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_point(data = top10_phylum.scaling1, color = 'blue') +
geom_text_repel(data = top10_phylum.scaling1, aes(label = group), size = 3, box.padding = unit(0.5, 'lines'), color = 'blue') +
labs(x = ca1, y = ca2, title = 'I 型标尺，仅展示 top10 丰度种群', color = NULL)

#ggsave('ca_phylum.pdf', p, width = 5.5, height = 5)
ggsave('ca_phylum.png', p, width = 5.5, height = 5)

评估有解读价值的CA轴

类似PCA，CA也只是一种探索性分析方法，主要目的是在低维空间尽可能多地展示数据的主要趋势特征。最终选择多少个观测轴、选择哪些轴并没有统一的标准，自己看着来。

通常的做法就是选取前2-4个排序轴，然后提取对象或变量在这些轴中的坐标，绘制排序图查看对象或变量沿这些轴的分布是否出现了某种规律。借此评估数据的趋势特征。

排位偏后的轴（通常就是第5轴之后）一般就不考虑了，因为它们的特征值太低。

可以通过柱形图由高到低展示所有的特征值，然后根据特征值做个评判，多少个排序轴值得保留和解读。（评估方法和PCA类似，也只是一种参考，并不是强制的）

#帮助评估有解读价值的 CA 轴
#Kaiser-Guttman 准则
#计算特征值的平均值，可据此保留高于平均值的主成分轴
ca_eig[ca_eig > mean(ca_eig)]

#断棍模型
#将单位长度的“棍子”随机分成与 CA 轴数一样多的几段，由高往低排序
#可据此选取特征根大于所对应断棍长度的轴，或者总和大于所对应断棍长度总和的前几轴
n <- length(ca_eig)
bsm <- data.frame(j = seq(1:n), p = 0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n+1-i))
bsm$p <- 100*bsm$p/n
bsm

#绘制每轴的特征根和方差百分比
par(mfrow = c(2, 1))
barplot(ca_eig, main = '特征根', col = 'bisque', las = 2)
abline(h = mean(ca_eig), col = 'red')
legend('topright', '平均特征根', lwd = 1, col = 2, bty = 'n')
barplot(t(cbind(100 * ca_eig/sum(ca_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('bisque', 2), las = 2)
legend('topright', c('% 特征根', '断棍模型'), pch = 15, col = c('bisque', 2), bty = 'n')

友情链接